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This is the fourth paper in a series that attempt to put forward a consistent framework for 
modelling solid regions in neutron stars. Here we turn our attention to axial perturbations of 
spherically symmetric spacetimes using a gauge invariant approach due to one of us. Using the 
formalism developed in the first paper in the series it turns out that the matter perturbations 
are neatly expressible in terms of a "metric" tensor field depending only on the speeds of 
shear wave propagation along the principal directions in the solid. The results are applicable 
to a wide class of elastic materials and does not assume material isotropy nor quasi-Hookean 
behaviour. The perturbation equations are then specialised to a static background and are 
given by two coupled wave equations. Our formalism is thus slightly simpler than the previ- 
ously existing results of Shumaker & Thorne, where an additional initial value equation needs 
to be solved. The simplification is mainly due to the gauge invariance of our approach and 
shows up also in somewhat simpler boundary conditions. We also give a first order formula- 
tion suitable for numerical integration of the quasi-normal mode problem of a neutron star. 
The relations between the gauge independent variables and the, in general, gauge dependent 
perturbed metric and strain tensor are explicitly given. 



1 Introduction 

1.1 Motivation 

In recent years it has become increasingly clear that, in order to extract information from observations 
of neutron stars, we must first understand the dynamical behaviour of matter beyond the perfect fluid 
approximation. Much effort have been put in to understand e.g. the dynamics of superfiuids including 
various viscous phenomena, see [T| for a recent review. This topic is essential in order to realistically 
discuss e.g. glitches and gravitational waves arising from secular instabilities. These events may also require 
understanding of magnetic interactions if the field is strong enough to influence the dynamics in the region 
where the matter currents are strong. This is certainly the case at least in the outer layers of magnetars 
and the outermost crust in pulsars but may even be important in the interior regions if the magnetic field 
is buried during the formation of the star. Glitches will require a thorough understanding of the behaviour 
of the rigid parts of the star and its interactions with normal and super-fluids. Recently a formalism was 
presented [1] suited for handling elastic solids permeated by a superfiuids, including, in a MHD like manner, 
magnetic fields. 

*E-mail: max@physics.muni.cz 
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The objective of this series of papers is to develop a coherent framework for purely elastic components 
in compact objects. In this, the fourth paper, we consider axial perturbations of spherical stars. It is worth 
pointing out that axial oscillations in neutron stars may already have been observed in the aftermath of giant 
flares in soft gamma ray repeaters, where quasi periodic oscillations with frequencies from 18 Hz reaching 
into the kHz range are observed [T3l [301 1231 HI EI] . Although the enormous magnetic fields believed to be 
present in these objects, which are well modelled within the so-called magnetar model [10) . will couple any 
torsional oscillations to the core within about one oscillation period via Alfven waves [T21 [H5] a recent toy- 
model suggests [12] that the frequencies will nevertheless be close to the purely elastic ones. Therefore, even 
in this type of systems, it is important to understand the elastic oscillation spectrum. This spectrum has 
recently been studied, within the relativistic Cowling approximation, in [26]. The local magnetic effects (i.e. 
neglecting the global nature of the modes) of the torsional mode spectrum has been studied for a poloidal 
dipole field [25] (see also [22]). These calculations were recently extended to include purely axial global 
oscillations [25] . 

The Newtonian theory of axial perturbations of elastic bodies is well developed, see e.g. [21] ■ The general 
relativistic problem was pioneered by Schumaker and Thorne [27] (hereafter ST) who developed a detailed 
theory including many relevant limits. Our treatment is different from theirs in several important ways. 
First, we use the gauge invariant perturbation formalism due to one of us [14] . This formalism assumes that 
the perturbations are axisymmetric, but as discussed e.g. in [6], this implies no restriction for spherically 
symmetric backgrounds. At a practical level, the difference between our and ST's approach shows up in 
the fact that our final equations consist of just two coupled wave equations and thereby dispenses of the 
additional initial value equation, being traceable to the gauge choice, in ST. Moreover, the formalism applies 
to any spherically symmetric (static or not) background. More importantly, our treatment does not from 
the outset assume an isotropic background which makes it applicable to the exotic ( "nuclear pasta" ) phases 
believed to exist in the lower crust of the neutron star (e.g. [24] ) . Another minor difference is that our 
formalism in the main part of this work is valid for non-linear clastic equations of state. That is, we do not 
assume a quasi- Hookcan equation of state from the outset. Practically this should not be very important 
since the behaviour of perfect elasticity (which we do assume) is likely to break down roughly when any 
nonlinear corrections to quasi-Hookean behaviour become important. 

As mentioned above, our treatment will not be limited to isotropic backgrounds. There are two funda- 
mentally different cases in which non-isotropy occur, either the material is intrinsically non-isotropic (as in 
the pasta phases) or the star is already strained. One may think that the latter is fairly unimportant since, 
at any given time, the crust will be close to its unstrained state simply because it will break otherwise. How- 
ever, the same limitations that apply to the background applies to the perturbations so that the background 
strain may well be of the same order as (or larger than!) the oscillations. An example where this may be 
important is provided by the flare model of Duncan [3] where the initial crust quake sends seismic waves in 
the (already strained) crust thereby causing it to fracture at other locations too. 

The treatment relies heavily on Paper I in the series [15] (hereafter Paper I), where the theoretical foun- 
dations of relativistic elasticity were outlined following the framework of Carter and QuintanaEj. For con- 
venience and in order to introduce notation a brief review of the formalism of Paper I is given below. 

1.2 Relasticity 

Any description of continuous media is, at least implicitly, based on the use of an abstract base manifold, the 
matter space, X say. This space (which in our case is three dimensional and Riemannian) can be thought of 
as a book-keeping tool which assigns unique labels to each fluid clement via a map that takes each flow-line 
on spacetime (M say) to a point in matter space. We may use this map to push-forward the contravariant 
metric g ab on M to A, 

9 AB = **S a& , (1) 

where abstract index notation E2] with capital letters on X and lower case letters on M is adopted. We then 
define the tensor r\~ 1AB to be the value of g AB that minimises the energy density at a fixed number density 
n. We now define tjab through the relation rf 1AC t\cb — S A b- Pulling back t\ab to M (rj a b = *S?*r]AB) we 



2 



may define the constant volume strain tensor according to 

Sab = -^{hab - Vab), (2) 

where h a b is defined by h ab = g a b + u a u b with u a being the unit tangent vector to the flow-lines. Thus, in 
effect, s a b measures the difference in geometry between the natural, unsheared state and the actual physical 
state. 

In this paper we will use the simplifying assumption that the elastic structure deforms conformally under 
pure compression. As discussed in Paper I this allows us to fully describe this structure by a fixed (i.e. 
n-independent) metric tensor field ftAB — n 2 ^rjAB on X. One may think of this metric as measuring the 
relative positions of the particles in a locally relaxed state. 

Using the fact that the eigenvectors of k ab = ^*kAB are eigenvectors also of the stress tensor p a b a 
preferred tetrad completed by the matter four-velocity u a may be introduced. Denoting the eigenvectors by 
e° we use Greek indices to numerate the space-like basis vectors. The eigenvalues of k a b correspond, in 
a loose sense, to (squared) linear particle densities, whereas the eigenvalues of the strain tensor are the 
pressures as measured by a comoving observer in the direction of e° . 

1.3 Perturbation theory and notational remarks 

We shall assume that there exist a family of solutions to Einstein's equations, parametrised by A sajQ We 
assume that the solution for a specific A (equal to zero say) to be known and expand around this solution. 
We take 6 in front of any tensor field to represent the value of d/dA on that field evaluated at A = 0. We 
shall have occasion to define perturbed quantities without the explicit perturbation symbol 6. The indices on 
such quantities are raised and lowered by the unperturbed spacetime metric, whereas care must be exercised 
when fields have retained the S, e.g. 

5e« = ^ 9 ab ^ b . (3) 

In order not to have too cluttered formulae we shall not notationally distinguish background fields from the 
full family fields but instead point out the validity of the equations in the text when confusion may arise. 

In order to conform with the notation in the preceding papers in this series [13 I17L I16j as well as with 
the work of Karlovini[14] on which much of this paper is based, we use some non-standard notation. In 
particular, the perturbed metric will be denoted by ^ ab and the flow-line orthogonal piece of the metric is 
denoted by h a b, i.e. precisely the opposite of the definitions in e.g. jSJ- 

We will use geometric units such that c = 1, G = 1 and the Einstein equations take the form 

Z ab := Rab - K{T ab - \T c c g ab ) = (4) 

where R a b is the Ricci tensor, T a b is the energy momentum tensor and k is the coupling constant which, in 
conventional units, takes the value k = 8irGc~ . 

2 General perturbations 

Perturbation theory of elastic media has already been considered by Carter in some detail in [5] . Since it seems 
computationally advantageous to use the eigenvalue formulation in many practical situations we shall devote 
this section to deriving the perturbations of the eigenvectors and from that deducing the perturbed stress 
energy tensor. The derivation is performed in an arbitrary (identification) gauge and the final expression 
reduces to the formulae of Carter when either a Lagrangian or a Eulerian gauge (given by some displacement 
vector field £ a ) is chosen. Since we are assuming a perfectly elastic conformally deforming material (so we 

1 We leave aside the technical, but important, issue of existence of such a family. On physical grounds, for the applications 
in mind in this paper, it is apparent that such a family should exist. 
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have a fixed matter space metric Uab) we shall consider the perturbed matter space and spacetime metrics 
to be our fundamental variables. 

The perturbed eigenvectors are easily derived from the identities 



g ab u a u b = -1, g ab u a e fl b = 0, gabe^e^ = 6^, 
k ab u a u b = 0, k ab u a e^ b — 0, k ab e^ a e u = n^S^^, 
which, when perturbed, yield 

u a u b j ab + 2u a 5u a = (5) 

u a e ll b lab + e M „ 8u a + u a <5e/ = (6) 

e /1 a e 1/ 7 6 + e Ma 5e v a + e va be^ = (7) 

u a u b Sk ab = (8) 

u a e^ b 5k ab + n 2 e^ a 5u a = (9) 

e^e^Skab +n^e fia Se u a + n^e va 5e^ a =6(^)6^, (10) 

where jab = Sg ab - Solving for Su a , <5e M a and 5(n?) gives 

3 

5u a = \u b u c lbc u a -Y, V 2 u V^fcc e M a (11) 
Se^ = — ^u b efl c (5k bc - n*-y ab ) u a - \ e^e^i^ e M a + ^ 2 - 2 e^ b e v c (6k bc - n^ bc ) e u a (12) 

= e^ a e, b (5k ab - n 2 lab ). (13) 
The perturbed basis one-forms are similarly found to be, 

3 

8u a = - \u b u c ~i bc u a - Y n^ 2 u b e c ^(Sk bc - n 2 ^ bc )e^ a (14) 

Se^a = -n~ 2 u b e c p: 6k bc u a + §e£e°7& c e M0 - ^ 9 ^ 9 e^e^(<5fc fce - n^7b c )e !/a (15) 

One may be worried that the difference between the linear number densities that appear in the denominator in 
the expressions (IT21 and (TT5"]) will cause divergences when there are degenerate eigenvalues on the background. 
However, as we will see, for the physical quantities relevant for this work, the sums are always convergent. 
These relations makes it a simple (but tedious) task to write down the perturbations of any tensor field in 
terms of the perturbed metrics. In particular the general perturbation of the stress energy tensor takes the 
form 

3 3 

ST ab = (5p - pu c u d ~f cd )u a u b + sXSpp + p (U e^e^7 C( i)e /10 e /1 6 - 2 ^ u c e^[n~ 2 (p + V^)^k cd - P7cd]w( a e jui) ) 

^i— 1 M— 1 

3 3 

+ e M e ^ [ n v 2 (P + P») v l±v( Sk c<i ~ n llcd) + Puled] e^ a e vV) (16) 

^i— 1 I^/i 

where v^v is the speed of shear waves along the /^-direction with polarisation vector along the ^-direction, 

{p + p v )v ±v = — 17 
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as derived in Paper I. It is now apparent that the isotropic limit is well defined. The perturbed energy density 
and eigenpressures can also be expressed in terms of the perturbed metrics by the use of (|13[) (remembering 
that the equation of state is considered to be given as a function of the linear particle densities) as 

s p = 2J2(p + p^ n » 2 ( Ska <>- n h"bK4 ( 18 ) 

5p » = \ X! n v^- n u 2 { 5k ab - nl~f ab )ele b u (19) 

(20) 

However, these will be seen to vanish for the axial case which is the main concern of this paper. One can 
confirm that, by comparing these relations to the expression for the elasticity tensor given in Paper I, one 
obtains the expression given by Carter [5], a calculation that is best done in a Lagrangian gauge for which 
u a Sk ab vanish. 



3 Axial perturbations 

We shall now specialise our considerations to axial (in a sense to be defined) perturbations. We use the 
approach of Karlovini 14 . For convenience, and to introduce notation, we briefly review this formalism 
below and then proceed to specialise it to perfectly elastic matter. 

Consider a spacetime with a Killing vector field r/ a . If this Killing vector field is nearly hyper-surface 
orthogonal, i.e. its twist 

n a = ^abcd^^ (21) 

can be considered to be small, we may treat the deviation from hyper-surface orthogonality as a perturbation. 
Note that this general considerations allows, as examples, studies of axisymmetric almost spherical systems 
as well as stationary almost static systems (with no required symmetry on the space-like hyper-surfaces). 
The main focus in this paper is the spherical case, so that rf is the axisymmetry generator and reduces to 
one of the SO (3) generators on the background. In order to make progress we split the metric according to 

g ab = ± ab +Fjj, a fi b , (22) 

where 

F = r) a r} a , f i a = F- 1 Va , r) a ± ab =0 (23) 

On a spherical background the function squared norm of the Killing vector has the form F = (r sm8) 2 . The 
polar and axial perturbations of the spacetime metrical are given by 

+ 7afc = {± a C U d + HaUbV^hcd = S± ab + (6F)fl a fl,b (24) 

~7 afc = 2T] c n {a ± b) d j cd = 2ri {a Sfi b} . (25) 

Taking the appropriate projections of Einstein's equations, 

± ac ± bd Z ab = (26) 
t/V Z ab = (27) 
± ab V c Z ab = (28) 

and using a partial identification gauge fixing amounting to 

Srf = (29) 



Note that the terms 'axial' and 'polar' will be inappropriate for perturbations of non-spherical spacetimes. 
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it is possible to show [T3] that, when evaluated at £l a — 0, the axial part of the metric (|25[) satisfies (|2U)) and 
(f2"T|) identically, whereas the polar part (|2"4"|) satisfies ((25)) identically. Hence these perturbations decouple 
and may be treated separately. 

We henceforth restrict ourselves to the axial case whence 

S± ab = 0, SF = 0, (30) 

so that 

lab = ~lab = 2V(a$Vb) (31) 

We next introduce the notation 

Q ab = 2V [a 6,i b] , J a = 26(± ah r) c T bc ) (32) 

as well as their (restricted) duals 

Q a = -±Fe abcd Vb Q cd , J ab = e abcd J c V d (33) 

Note that both the perturbation vectors Q a and J a (as well as their duals) are gauge invariant since the 

corresponding background quantities vanish. In terms of these fields the full perturbation equations takes 
the Maxwell-like form [T3] 

V b (FQ ab ) = K.J a , (34) 

V a J a = 0, (35) 

or in the dual form 

2V [a <2 6] = -nJ ab , (36) 

V a (F- 2 Q°) = 0, (37) 

V [a J bc] = 0. (38) 

To proceed it is evident that we must compute the axial matter current vector J a . To this end we assume 
that one of the eigenvectors of the matter space metric, eg say, is aligned with the axisymmetry generator 
on the background, i.e. 

e 3 a = F- 1/2 r) a . (39) 

It is worth pointing out that we hereby restrict ourselves to the case when the Killing vector is space-like, 
but that no loss of generality is implied in the spherically symmetric background case. Up to this point the 
discussion has been valid (with appropriate changes of terminology) to background spacctimes admitting 
any hyper-surface orthogonal Killing vector. 

To find J° we project (fT6|) with ±_ ab rf and use the relation 



c i db ac bd ( a n\ 

o^- = -.9 9 led, (40) 

so that, after some algebra, 

r 1 c T bc S± ab = -p 3 r ] c ± ab i bc , (41) 

-L ah V c 5T bc = nf(p + p 3 ) (^-u a u b + ^ e »U^j ^ ~ n hbc) + vrf ± al 'ibe (42) 

= (p + P3 )FS ab K b ~ f] c T bc 6± ab , (43) 



where S is the "metric" 



S ab = -u a u b + J2v,\ 3 e, a e, b , (44) 
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and K a is the gauge invariant (vanishing when unperturbed) quantity 

K a = (Fn£)- 1 6(± a b r l c k bc ) = (Fn 2 )- 1 _L a b r) c (5k bc - n 2 lbc ) = {Fn 2 )- 1 ± a b V c Sk bc - Sp a . (45) 
Combining terms we finally find, 

J a = 2(p + p 3 )FS ab K b . (46) 

The tensor S ab is vaguely analogous to the so called acoustic metric (sometimes denoted G ab ) used in e.g. 
the study of analogues of black holes in fluid mechanics in the sense that it is related to the propagation of 
waves (although, in this case, it is shear waves rather than sound waves). 

We still need to evaluate K a . We therefore make the very natural choice that the matter space is 
axisymmetricH and introduce material space coordinates (x 2 , </>), i = 1,2, such that 

V a V a x l = 0, ifV a i = 1. (47) 

We also assume that the pull-back of the material space metric takes the form 

k ab = KaVa&VtfP + FV a 4>Vb4>, (48) 

hi 

where the metric components and F depend on x l only. Contracting this relation twice with e 3 we find 
that the function F may be expressed as 

F = Fe%\e\k ab = Fn 2 , (49) 

which holds on the background. This allows us to evaluate K a to be 

K* = V*8<j>-6na. (50) 

We may note in passing that an equivalent form of the perturbation equations may be found directly in 
terms of the gauge invariant one- form K a . The equations of motion are then 

V b (FV^ b K^) = K(p + p 3 )FS ab K b , (51) 

together with V a J a = with J a given in terms of K a by (|4l)|) . We shall, however, not pursue this form 
further. 

3.1 Spherically symmetric background 

We shall now specialise the equations further by assuming that the background is spherically symmetric. To 
ease the presentation we start by introducing some notation. As much as possible we shall use the notation 
of Paper II [17] in this series. The background metric will be decomposed as, 

9ah = jah + t a b, (52) 

where t ab is the metric on the hyper-surfaces spanned by the SO(3)-generators and is taken to be represented 
by the line element 

ds 2 t =r 2 (d0 2 +sin 2 6»d</» 2 ), (53) 

where r is the Schwarzschild radius. Introducing the shorthand notation ? ,£l = e", the orthogonal two 
dimensional Lorentzian metric is given by 

jab = -u a u b + r a r b , (54) 



3 One could, in principle, consider a non-axisymmetric matter space and constrain the mapping in such a way that spacetime 
is still axisymmetric. This feels highly unnatural, however, and we will not consider it further. 
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and the associated volume form is 

tab = -2u [a r b] (55) 

The covariant derivative operator associated with j ao will be denoted by T> a . Since the angular background 
eigenvectors will now correspond to degenerate eigenvalues we may further introduce the notation for the 
basis indices, 

p = 1 — ► fj, = r, ;ie{2,3}->^ = t, (56) 

that is, for example 

Pl=Pr, P2=P3=Pt- (57) 

We now proceed to separate out the angular dependence from our perturbations equations. The discussion 
will be very brief and the reader is again referred to [H] for details. 

Using the dual form of the perturbation equations it is evident from the closedness of J ao expressed in 
that we may introduce the vector potential Y a according to 



J„i = 2V [o y 6] . (58) 

Furthermore, this vector can be taken to be orthogonal to the SO(3)-generators. Proceeding to equation 
([3"Bf it is seen that 

Q a = V a $ - nY a (59) 

for some axisymmetric scalar <£>. Separation is achieved by putting 

$ = C(9)ri>, Y a = C(9)e ab J b . (60) 

The vector J a is a two-dimensional object related to J a by 

J a = (r 2 sin0)- 1 [C7 / (0).7 a - C{e)J(d/d6) a ], (61) 

for the two-scalar J constrained by ([55]) to take the value J = T> a J a . It is convenient to introduce the 
invariantly defined mass function m through 

9m 

1 = V a rV a r, (62) 

r 

as well as the function r invariantly defined on spherically symmetric spacetimes by 

R a b if = i K T?f, (63) 

i.e. hlvr is the eigenvalue of the Ricci tensor with respect to the eigenvector rf. A clearer physical inter- 
pretation is seen from the background Einstein equations which implies that r is the minus the trace of the 
2x2 block of the energy momentum tensor orthogonal to the SO(3)-orbits. In the present case that is 

r = p-pr (64) 

Noting also the relation 

2m 1 

rV a V a r = -nr 2 T (65) 

r 2 

the perturbation equations can straight forwardly be confirmed to reduce to the two-dimensional equations 

(V a V a - U)ip = kS (66) 

V a J a = J, (67) 

where 

S = re ab V a (r- 2 J b ), (68) 
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and the potential U is given by 

^=2 KT ~ — + ( 69 ) 

where Z is the usual harmonic index coming from the separation of variables. The angular equation is solved 
by setting 

C(0) = G;V 2 (cos 6) (70) 

where G l H 2 (y) is an ultra-spherical (or Gegenbauer) polynomial. 

In our case J a is given in terms of background fields and K a , where, as implied by eq. |5U|) . 

2V [a K b] = -2V [a <5/i 6] = -Q ab . (71) 

This means that Sfi a must be of the form 

C{e) 

(Z + 2)(Z-l)r 2 sin a 

Sfi a = r-\ a b Q b + V a / G {I = 1) (73) 

where 

Q a =V a {r^)-Ke ab J b , (74) 

and /g is a free function (which we will see in section [5] correspond to the choice of gauge) . This in turn 
implies that K a must be of the form 



S ^a = - 71 —, \ 5 . 3 £a b Qb + V a / G (i > 2) (72) 



X- = ^e a b Q b + 1 7 V a fsin- 3 6»C"(6»)r-Vl 

(Z + 2)(Z-l)r 2 sin 3 # (Z + 2)(l-l) 1 ' J 



= n^n' { °l 2 ■ 3 J ta b Q b + r 2 V a (r-^)} -^L-VVa*. (/ > 2) (75) 
(t + 2)(t — ljH sm sin 6* 

^a = r- 2 e a b Q h + I? a (r-V) (i = 1) (76) 

for some spherically symmetric scalar Setting 

5 o6 = -u'li 6 + t&rV, (77) 

we have 

S ab = S ab +v t l t e 2 a e b . (78) 

For / > 2 we then find 

J Q = 2(p + Pi )^ ah ^ 

2(p + ft ) \ C'(9) ^[^Q^^^i^^^^-i^^J (7g) 



sin (9 [(l + 2)(l-l) 
Comparing with eq. ()61|) . we obtain 

[f 5 a 6 + (i + 2)(I - lb'M^J 6 = £ 5 ah [e 6 c P c (r^) + r 2 2? b (r- V)] (80) 

«J = f « t i t r~V (81) 

where 

£ = 2«r 2 (p + Pt )- (82) 

Solving for nj a we find 

K J« = (~C u a u b + C ir a r b )[e b c V c (r^) + r 2 V b (r- x V )\ (83) 
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where 



C = 
Ci = 



£ 

£+(l + 2)(l-l) 
£v? ± + (l + 2)(l-l)- 



When I = 1 a similar analysis shows that J = and that J a is given by the expression 



with 



Co = 


£ 


£ + 1 


C x = 




£v r \ + l 



Finally, rewriting the perturbation equations (|r]rj)) - ((^7)) in terms of the quantity Q a defined in ([74")) 



find 



84 
85 

8G 
87 
we 

88 
89 



where Q a , J a and J are given in terms of the two-dimensional scalar fields ip and <p in eqs. (fT4"|) , ((55)) and 
(|5T|) respectively. These equations are valid for any I > 1 provided one remember to put J = when i = 1. 
The non-radiative nature of the / = 1 case is hinted by the form of the equations which are just decoupled 
conservation equations for the two "currents" r~ 2 Q a and J a . 

In order to better understand the content of these equations it is useful to write them out in a suitable 
coordinate system. For brevity we will from now on only consider the radiative case I > 2. It is advantageous 
to use explicitly conformally flat coordinates given by 



u a = -e u (&t) a , r a = e iy (dr*) a . 



(90) 



which will reduce to the usual Regge- Wheeler radial gauge on static backgrounds. Next we introduce the 
notation 



for any X a . Then we have 



X t =u a X a , X r =r a X a 

Kj t = e- v C [(rip)' +r 2 (r- 1 cpY 
Kj r = e- u d UripY +r 2 (r- 1 ip)' 

Qt = e-" [(1 - CiXnfO' - cV(r-V)' 

Qr = e- v [(1 - C ){ri>)' - C r 2 (r-V)' 



where dots and primes denotes derivatives with respect to t and r* respectiveljQ- Using 

W t =r~ 1 e v Q t , W r =r- 1 e»Q r 



(91) 

(92) 
(93) 
(94) 
(95) 
(96) 

(97) 



'This notation should not be confused with the use of dots and primes in Paper II |17| . 
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as auxiliary variables we can cast the perturbation equations as a first order system, 



f r' ,, L 



-Wt + W r + -Wt--Wr-e 2u ^ = S (98) 
-W r + Wi--W r + -W t + e 2l/ ^±(p = (99) 

-L(rip)' + f w r 2 x r 2 (r-V)' + (£u r x + i)rW f = (100) 

L(r^)' - £r 2 (r~V)' ~ (£ + £)rW r = (101) 

where L = (I + 2)(Z — 1) and we have used the definitions of Co ([54")) and Ci (|55|) . We may note here that 
the principal part of these equations decouple into two systems, one in ip and tp which has a characteristic 
propagation speed v 2 j_ and one in W t and W r whose speed is 1. We may naturally interpret this as the 
existence of two families of modes, the shear modes and the axial gravitational ui-modes. 

This is as far as we go in the general case. Note that this system is suitable for numerical integration since 
it does not contain derivatives of background quantities that can be discontinuous. As we will see below, it 
also makes maximal use of the junction conditions so that W t , W r and ip are all everywhere continuous. 



4 Static background 



The above derived system (|98[) - (|101|) of equations is valid for any spherically symmetric background and 
could therefore be applied to e.g. collapsing bodies. However since our interest here lies on matter with 
non-zero shear modulus such configurations are not realistic (since the strain would inevitably grow beyond 
the breaking strain of the material). One could of course study axial perturbations of a radially oscillating 
star, but henceforth we will restrict our attention to the case where the background is static. We can then 
drop all dotted background quantities. It is then possible to combine - (jlOOj) to give a wave equation 
for Wt sourced by (p. In addition the equations ([99jl - (|10ip can be combined into a wave equation for ip 
with Wt as a source; 



£v 2 ± 1(1 + 1) 



e^£(v r \-v t \ t ) 



-1 / 









'■p 



-p + 



£r r 2 



cr £ 



(102) 
(103) 



Of course, many possible reformulations exist of this system of equations if new dependent variables are 
chosen. One such reformulation would have been to choose tp instead of Wt as the first dependent variable. 
This would have followed the general guidelines in [13] more closely, but would not have lead to simplified 
formulae in this case. We do not feel that any possible second order system reformulation (that we have 
found) is superior in any important ways to (|102[) - (|103p . In any case, for numerical integration it is preferable 
to use the first order system (p?8")) - (|101[) . We could also give more explicit expressions of the terms involving 
derivatives of background quantities by making use of the results of Paper I. However, either one then keeps 
the elastic equation of state free in which case third order partial derivatives of the energy per particle will 
appear, or one specify e.g. a quasi-Hookean equation of state leading to rather lengthy expressions. We 
do not feel that this leads to any further understanding of the problem at hand. Instead we give explicit 
expressions below suitable for numerical integration. We shall therefore be satisfied with this form of the 
wave equations and pause a moment to discuss some of their properties. 

First, it is clear that for vacuum (|103[) is trivial and (|102[) is just the usual Regge- Wheeler equation. For 
this reason we may refer to (|102j) as the gravitational wave equation. Second, if we take the isotropic limit 
we have v 2 ± — vf ±t = 2nr 2 fi£~ 1 1 where jj, is the shear modulus, so that the source terms in (|102|) reduce to 



2 K (e 2l/ /i)V 



(104) 
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so that we see that, as pointed out by Dyson for weak gravity, e" ~ 1 ~ constant, gravitational waves 
couple only to the gradient of the shear modulus. Another important case is the perfect fluid limit. Here we 
see that the system reduces to 

-W t + W t " + e 2 " - \n{p -p)- fcil) Wt = (105) 

-r(p" = (rWt)' (106) 

The first of these equations is the standard result, see e.g. [3 El HI]. The second equation is now redundant 
since the variable ip does not encode any relevant physical information in this case[l. It is in fact easy to 
show, using ([92 ]) -([94 | and (|106p . that the two dimensional matter current J a satisfy J r = 0, Jt = so that 
the well known result that axial perturbations of perfect fluids reduce to stationary currents is evident. 

4.1 Boundary conditions 

The perturbation equations must be accompanied by suitable boundary conditions. These conditions were 
treated in detail by Karlovini for general matter sources |14j and we merely state the results here; At any 
boundary (including e.g. interfaces between different layers or the surface) we require that the one-form 
Q a as well as the scalar ip should be continuous. These conditions follow from the requirement that the 
first and second fundamental form on any hyper-surface with normal n a (say) are continuous. In particular 
this implies that the traction n a T a i, is continuous as well |20j . In ST the condition of continuous traction 
was enforced as an additional constraint. In view of the above it is clear that this is only necessary as a 
gauge condition relating the gauges on either side of the hyper-surface. An additional boundary condition is 
provided by requiring regularity at the centre. Finally, if we consider an isolated object, we need to impose a 
condition of outgoing gravitational waves at infinity. We discuss these boundary conditions in more detail in 
the next section where we give a computational recipe for quasi normal mode solutions to the perturbation 
equations. 



5 Computational algorithm 

We shall here present a procedure for finding quasi-normal mode solutions to the perturbation equations in 
situations relevant for neutron stars with a fluid core, a solid crust and, possibly, a fluid ocean. We therefore 
assume harmonic time dependence, i. e. that the temporal dependence of the independent variables is given 
by e ZUJt . We have found it slightly preferable to present the equations in Schwarzschild coordinates, related 
to the Regge- Wheeler coordinates through 

e"dr* = e A dr (107) 

where A is given by 

er 2X = 1 108 

r 

When we wish to solve our perturbation equations numerically it is preferable that all variables scale with 
the same power of r near the centre. To obtain this behaviour (as well as some other features) we redefine 
the variables in the system (j9"B")) - (|10ip (specialised to a static background) according to 

X x = r-^Wt (109) 

X 2 = -iLu r- l W r (110) 

X 3 = -iw e ! V-'- 1 V> (111) 

X A = u>le v r- l <p (112) 



5 The extra equation correspond in a sense to the freedom of choice of mapping between M and X. For perfect fluids the 
labelling of fluid elements plays no role in the dynamics so there is some degree of degeneracy in defining the map. 
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where luq = e~ Uc uj and v c is the central value of v. The reason for introducing these variables is that if one 
wish to solve the background equations simultaneously with the perturbation equations using a shooting 
algorithm the central value of v is unknown since it is determined by matching the background solution 
to the exterior Schwarzschild solution at the surface of the star. The above scaling allows for writing the 
equations in terms of loq and v§ = v — v c and, once the solution is found scale it back to physical units. 

It is easy to show that for perfect fluids the last two of the first order equations become algebraic relations 
amounting to 

X z = -e v °X x X 4 = -e u °X 2 (113) 

whereas the other two become just 



r— ± = -(I + 2)X X 
dr 



z x - Ua X 2 



„-2Z-l 



dr 

Disregarding the singular solution scaling as 
centre according to 

Xi = Xi+Oir 2 ) 

X 2 = -{l + 2)X 1 +0{r 2 ) 

X 3 = X x +0(r 2 ) 

X 4 = -(l + 2)Xi +0{r 2 ) 

where X\ is an arbitrary constant. The equations in a solid becomes, 

,£v 2 



(114) 
(115) 



near r = it is easy to see that the variables leave the 



(116) 
(117) 
(118) 
(119) 



= _ (i + 2 ) Xl e x ^X 2 - e x ^f 
dr r l io 







£v r± r 



dX : 
dr 
dX 3 
dr 
dX 4 



2 A— i/o 2 2 

— — g . ■ i ■ 



t4r*Xi -{I- l)X 2 + e LX$ 



= e x ^-^X 2 + [rv, r -(1 + 2)}X 3 + e A ^ yX 4 

Li 1j 



dr 



-eVV 



where 



%[£v? x + L)Xx - e x - V0 r 2 uj 2 LX 3 + £v? x [rv, r -(I- 1)]X 4 
dv m+ T;nr 3 p r 



(120) 
(121) 
(122) 
(123) 

(124) 



dr r(r — 2m) 

and the shear wave speeds depend on the equation of state according to (fT"7|) . For the quasi-Hookean equation 
of state discussed in Paper I they take the form 

c 2 - 2 

£v rX = K/ir 



£v 



t±t 



2n(ir 



r^2 



"7 
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In vacuum the equations are 



X 3 = -e»°X u 



dX x 
dr 
dX 2 



-(l + 2)X 1 



L — 



X 4 = unconstrained 

re" c 



r - 2M 
? 2l/c oj 2 r 3 
r - 2M 



-X 2 



Xi — (I — 1)X 2 



(125) 
(126) 

(127) 

(128) 
(129) 
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where M is the total mass of the star. 

Using the expansion (|116[) - (|119p we may now integrate the fluid equations to the crust core interface, 
where we need to impose the boundary conditions. In the present variables these are just the continuity of 
Xi, X2 and X3. Note that X4 is free at the boundary and can in principle be set to any value. However, as 
we proceed with the integration in the solid we eventually encounter the next boundary. If this boundary is 
an interface to a fluid phase or vacuum we have two conditions on the variable X% , i) it must be continuous 
and ii) it has to satisfy the constraint (|113p (which is identical to (I127[) in the vacuum case). In general 
these will not be satisfied signalling the wrong choice of value of X4 at the previous boundary. The remedy 
to this situation is to find a second independent solution to the equation in the solid phase and then to take 
the linear combination of the two solutions which satisfy all the boundary conditions. The simplest way 
to find the second solution is to start with X\ = X-x = X3 = and X4 — any value. Clearly, adding any 
multiple of this solution to the original one will not violate the initial boundary conditions so it is a simple 
task of solving an algebraic equation to find the correct total solution in this part of the star. It should be 
evident how to work out the entire solution to the problem using this algorithm regardless of the number of 
interfaces in the star. Finally, the quasi-normal mode frequencies u) are determined by making sure that the 
vacuum solution is described by outgoing gravitational waves. This specifies the solutions up to a scale. 



6 Identifying the metric and the strain 

Since neither the metric nor the strain on a strained background is gauge invariant we need to specify the 
gauge in order to evaluate them. We have already used some of our gauge freedom by setting 5if — 0, so 
any further gauge transformation has to be generated by a vector field, £ a say, that is Lie-dragged by the 
axisymmetry generator, 

c v c = -c c v a = (130) 

We may decompose C° according to 

<T = CI + iW\ ClVa=0, V a VafG=0 (131) 



Now, since 



C af x a = (132) 



due to the fact that i] a is a hyper-surface orthogonal Killing vector on the background it is clear that gauge 
transformations generated by a vector field orthogonal to, and Lie-dragged by rf does not affect the axial 
perturbations. Likewise, it is easy to show that a gauge transformation of the type ( a — fGV a does not affect 
the polar perturbations. Hence we need only consider gauge transformations of the latter type. Under such 
a transformation we have 

S^a -> 5[i a + V /g, (133) 
84>^5$ + f G . (134) 

Since r/ a \7 a S4> — Owe may use this freedom to set 5(f> = (fj This gauge choice is of a comoving, or Lagrangian 
type and is the natural way to measure the strain since e.g. the breaking strain will invariably be calculated 
in a comoving frame. We thus find that in a Lagrangian gauge we have A0 = so that /g = —5<p and 
A/z a = — K a , where we use A to indicate that a special gauge has been chosen. Using these relations we find 

As afc = -T] (a K b) + u {a Au b) = -f] (a h b) c K c . (135) 

The perturbed metric is just 

Ag ab = -2 V{a K b) (136) 



6 Note that Schumaker and Thorne uses the gauge freedom to set the 0-component of Sfi a = 0. This would correspond to 
setting /q = constant. It is also worth pointing out that, since these authors consider isotropic backgrounds, their perturbed 
strain tensor is gauge invariant. 
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and 

Ak a b = 2k bc ^ c K a \ (137) 

Usually estimates of the breaking, or yield, strain is given as a dimensionless strain angle, 9 max say. 
The precise relation to the components of the strain tensor depends on the type of deformation and the 
microscopic structure of the solid material. However, for simple (e.g. isotropic or cubic) structures under 
simple deformations (e.g. pure twist or shear) the nonzero components of s ab have the form 

s ab « |e. (138) 

Given the uncertainty in the literature on the value of max |25j it does not feel meaningful to digress too 
deeply into the subject of breakdown of elasticity. It is clear that the approximation of perfect elasticity 
will break down before the material actually cracks, so for all purposes of this paper we may assume that 
something catastrophic happens when any component of s ab exceeds some value of the order of ^O max - 

An interesting observation is that, even though the first order perturbations of the linear particle densities 
are all zero, so that the strain scalar also is unperturbed to first order, we may still estimate the second 
order contribution to the strain scalar and thereby estimate the energy stored in the elastic material. We 
may namely compute the strain scalar by first evaluating the total 

k a b = °k a b + Ak a b + A 2 k a b (139) 

for some unspecified second order perturbation A 2 k a b . Noting that the non-orthogonal pieces will not alter 
the strain scalar (all such pieces are contracted with orthogonal tensors) and using bold face letters to denote 
the 3 x 3-matriccs k = h a c k c b we may write (see Paper I) 



s 2 = 



1 



36 det k 



[(Trk) 3 -Tr(k 3 )-24detk] (140) 



It then turns out (after some straightforward but tedious algebra) that the second order perturbation terms 
only comes in multiplied by the background strain. These terms will therefore always (due to the small 
numerical value of the breaking strain) be much smaller than other second order terms. Neglecting them 
leads to an expression of the form 

s 2 = si + X -Fh ah K a K b + ... (141) 

where so is the background value. 

In summary, once a solution to the perturbation equations are found using the procedure described 
in section [5] it is easy to find the metric and stress tensor perturbations from equations (1136|) and (1135p 
respectively. Since the solutions will only be defined up to a scale one may then set this scale such that the 
largest component of s ab + As ab is smaller that some maximum value of order ^O max to find the maximally 
allowed amplitude of the perturbations consistent with elastic response. It is also straight forward to estimate 
the maximal stored elastic energy density by the formula 

Aplastic = fis 2 (142) 
where jj, is the shear modulus and s 2 is given by l|14ip . 



7 Conclusion 

We have developed the general relativistic theory of torsional oscillations in elastic matter. In the light of 
the recent very exciting observations of quasi-periodic oscillations in the tail of giant flares in soft gamma- 
ray repeaters, with frequencies matching well the expected spectrum of such modes, we believe that it is 
important to have a well founded theory for these modes. We have improved on the previously existing theory 
(mainly ST) in several respects; The equations presented here are gauge invariant and are valid to second 
order in strain for any equation of state as long as the material is of a conformally deforming type. Various 
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gauge-issues have been resolved and resulted in simplified equations and boundary conditions compared to 
ST. The elastic response to perturbations is conveniently parametrised in terms of the shear wave velocities 
and should be straight forward to apply to the anisotropic "pasta" phases near the crust-core boundary of 
neutron stars once the relevant elastic parameters have been estimated. 

The results of this paper has in fact already been applied in the relativistic Cowling approximation |26| . 
Although it is expected that the spectrum should not change substantially when gravitational degrees of 
freedom are included it is nevertheless important to check this numerically and work in that direction is in 
progress. 
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